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TECHNICAL MEMORANDUM 


A GENERALIZED REUSABLE GUIDANCE ALGORITHM FOR 
OPTIMAL AEROBRAKING 

I. INTRODUCTION 


The use of planetary atmospheres for aerobraking as a means of saving propellant has been dis- 
cussed and studied extensively since the first paper on the subject appeared in 196 1. 1 Many studies have 
concluded that significant propellant savings can be realized 2-5 compared with all-propulsive orbit trans- 
fer. The amount of savings depends on the mission scenario and the weight of the thermal protection 
system (TPS) required to protect the spacecraft from the aerodynamic heating incurred during the deep- 
est penetration of the atmosphere. Reference 6 is an excellent survey of aerobraking studies. 

In the past, aerobraking guidance algorithms 7 have been developed with computational effi- 
ciency as one of the major design goals, an obvious carryover from an earlier period when computing 
capability was a mere fraction of what it is today. Thus, simplifying assumptions were made in order to 
avoid numerical integration of the equations of motion. As a result, the guidance algorithms were not 
very general or adaptable for future aerobraking mission planning. The speed of today’s computers 
permits more realistic, and therefore more general, real-time guidance modeling. A big advantage in this 
approach is the reduction in premission analyses required to determine appropriate values for mission- 
dependent parameters that are needed because of the simplifying assumptions. For example, the guid- 
ance algorithm for the aeroassist flight experiment (AFE) 8 used the following guidance law for deter- 
mining commanded bank angle: 


COS (<Pcmd) — COS (< p eq ) + K q 




( 1 ) 


where ^ is the bank angle required to maintain “equilibrium flight,” K q is a premission determined gain 
on dynamic pressure, .is a premission determined gain to provide guidance margin, and AT/ is a 
premission determined gain on altitude rate error (actual minus commanded). The gain, Kf, has two 
empirically derived values, one for the entry phase (also, the equilibrium glide phase) and one for the 
exit phase which begins when velocity magnitude has decreased to below a premission determined 
“trigger velocity.” Yet another premission determined “transition velocity range” dictates how to 
smoothly transition to the exit phase. There are also premission-determined values for the plane 
controller logic. The form for equation (1) was presumably derived on a highly empirical basis as are the 
gain values and some of the other mission-dependent parameters. The algorithm to be described in this 
report represents a more direct “mapping” from the real-world problem of successfully aerobraking into 
the aerobraking guidance domain where that problem is solved. The result, it is believed, is an algorithm 
that is less empirical and easier to use and understand. 


The next section will serve to introduce guidance concepts by discussing a simple apoapsis con- 
troller. Section III will discuss an optimal real-time apoapsis controller and the optimization results on 
which it is based. In section IV, a newly developed orbital plane controller will be discussed. Numerical 
testing and results are discussed and presented in sections V and VI. Finally, conclusions are made in 
section VII. 



n. BASIC REAL-TIME APOAPSIS CONTROLLER 


The strategy here is to numerically simulate, within the guidance, one or more planar aerobraking 
trajectories flown at constant bank angle. Based on the results of the simulated trajectories, the constant 
bank angle required to attain the target apoapsis radius at atmospheric exit can be determined. The 
following set of planar equations of motion, 9 for example, may be used for computing integrated 
trajectories: 


dr/dt = v sin y 
dv/dt = -Da- g sin y 

dy/dt = La cos 0/v + (v/r-gA>) cos y , (2) 

where r is vehicle position magnitude, v is inertial velocity magnitude, y is the inertial flight path angle, 
D a is the aerodynamic drag acceleration, La is the lift acceleration, and g is the gravitational acceleration. 
Estimated density as a function of altitude is needed in the integration process and can be computed 
from a simple onboard exponential model combined with a density estimator. Density, p est , is estimated 
each guidance cycle from navigation measurements 


Pest=2mD meas /(SC D V*) , (3) 

where Dmeas is the measured acceleration component in the relative velocity direction. From this, a mul- 
tiplier, krho is computed 


krho — Pest I PboddW » W 

and, thus for the current guidance cycle, the entire density profile can be estimated as 

p(h) = krho Pmodel(^) » 

the assumption being that the relative difference between the model density and encountered density is 
constant or at least slowly varying. In practice, to reduce the effect of noise, the krho used in equation (5) 
is actually the output of a low-pass filter driven by the k T ho in equation (4). Note that accurate calculation 
of D a and L a in equation (2) requires knowledge of relative velocity magnitude which is not obtainable 
from equation (2). Past experience shows that we can generally assume that the difference between 
inertial and relative velocity magnitudes is a constant throughout the trajectory. Alternatively, more 
generality could be gained by using a sixth-order differential system for the cost of additional 
complexity and computation. In any case, the problem of determining the appropriate bank angle, 0, is 
mathematically that of determining the zero of the nonlinear function 

^( 0 ) \r = rex “ rat-r a (r,v,y)\r = rex * ( 6 ) 

where rat is the target radius of apoapsis, and r a is the radius of apoapsis, a function of the simulated 
vehicle state, r,v,y. The initial conditions are derived from the onboard navigation system. There are a 
number of techniques available to solve equation (6) including Newton’s method or variants thereof. In 
practice, one Newton update (of required bank angle) is all that is required per guidance cycle, meaning 
either one trial trajectory plus one variational trajectory (in the case of the variational form 10 of 
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Newton’s method) or two trial trajectories (in the case of finite difference form of Newton’s method) 
must be generated per guidance cycle. (Actually, it is even possible to solve equation (3) using only one 
integrated trajectory per guidance cycle and information from the integrated trajectory from the previous 
guidance cycle.) The advantage of the variational form is that finite difference error does not affect the 
accuracy of the Jacobian used in the Newton update, whereas a disadvantage of the variational form is 
that anytime a change is made in the differential equation model, the variational equations have to be 
changed also — a tedious task. The disadvantage of the finite difference form is that robust logic has to be 
developed to generate accurate Jacobians, and that two separate trajectories have to be generated per 
guidance cycle. The advantage is that once the logic has been developed, it is valid for any guidance 
model changes that may be deemed necessary. 

To improve numerical conditioning, the solution of equation (3) is actually obtained by iterating 
on the cosine of 0 instead of 0 directly. If the solution for a given guidance cycle happens to result in 
Icos (0)1 > 1 (this is referred to as guidance saturation), then one simply limits the bank command to 0° 
or 180° as appropriate. For the next guidance cycle, one uses the solution from the previous guidance 
cycle as the initial guess even if Icos (0)1 > 1. TTiis is perfectly acceptable in computing integrated trajec- 
tories because, by inspection of equation (2), cos (0) is just a scalar multiplying die lift acceleration. 

A. Variational Equations 

A set of variational equations used to solve equation (6) is developed here. The Newton iteration 
for cosine of bank, c0, can be written 


c0 ,+ 1 = c0,- , (7) 

dyrfd(cQ) 

where yris the constraint function on apoapsis, and dy/fd(c<f>) is the sensitivity or Jacobian, both of 
which are evaluated at atmospheric exit. Using conservation of energy and angular momentum, the con- 
straint function can be expressed as 

yKnv.Y) = l+2/i(l/r a ,-l/r)/v 2 -r 2 cos 2 y/rj-. (8) 

The Jacobian, dy//d(c(t>), is expressed as 

dw dw dr dy dw 3v dy dw dy 

— — = — — — + — — + -f- — — . (9) 

d(c<j>) dr dyd(c<(>) 8v dyd(c<l>) “Y d(c<p) 

It is a routine matter to obtain the partials of y/ with respect to r, v, and y To obtain the rest of the com- 
ponents of the Jacobian, we make use of the relation 

d(dx(t)/da)/dt = d(dx(t)/dt)/da = d(f[x(t)))/da , (10) 

where dx/dt =f{x(t)). According to equation (10), the time derivative of the sensitivity of x (/) to some a 
is simply the partial of dx/dt with respect to a . Using this gives the following variational equations 
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d(drfdy)/dt = d(dr/dt)/dy= d(v sin y)!dy= dv/dysm y+v cos y , 
d(dv/dy)/dt = -D a (dp/dh dr/dy/p+2 dvfdyl V r )+2g dr/^yan ylr-gcos y , 

d(dy/dc$)/dt = LJv , 

with the initial conditions 


( 11 ) 


dr/dy(t = lb) = 0 

dv/dy(t = to) = 0 (12) 

dyfdcf (t = to) — 0 . 

Thus, we can obtain, via numerical integration, the partials, drldy, dvfdy, and dy/d(cp), at any given 
time. In practice, we simultaneously propagate the states, r, v, y, and the partials drldy, dy/dy, and 
dy/dcQ, to atmospheric exit resulting in all the information needed in the Newton iteration equation (7). 
Though it is not strictly correct to propagate the state and the partials simultaneously (the equations are 
to be evaluated along the nominal trajectory, equation (10)), doing so is a reasonable approximation and 
provides for efficient computation. The formulation and implementation of the variational equations can 
be numerically verified by comparing the Jacobian computed using finite differences and the Jacobian 
computed using the variational formulation. 

B. Numerical Integration 

Man y efficient numerical integrators exist today. The adaptive step size fifth-order Runge-Kutta- 
Fehlberg numerical integrator was chosen because it is accurate, efficient, and requires a relatively small 
amount of code and overhead. The integrated trajectories do not need to be computed very accurately 
because navigation errors and modeling errors cause loss of trajectory accuracy anyway. Experience has 
shown that 0.1 km in position error in the integration is adequate. Typically, no more than 15 to 20 
integration steps are required to compute an integrated trajectory. 

C. “Trappedness” Indicator 

To ensure that the integration process inside the guidance is robust, a test must be implemented 
to detect simulated vehicles that are “trapped” in the atmosphere due to too much negative lift being 
modeled. It is desirable to detect trapped guidance model vehicles as early in the integration as possible 
to avoid wasting computation time since trapped trajectories yield no useful information except for pos- 
sibly providing a lower (or upper) bound on required bank angle. Also, for the sake of robustness, the 
test must always detect trapped trajectories, and must never provide a false indication of “trappedness.” 
The following condition has been found to be a useful indicator of trappedness 

y<0 and dy/dt <0 . (13) 
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III. OPTIMAL REAL-TIME APOAPSIS CONTROLLER 


Extensive aerobraking trajectory optimization work has been done in recent years. 3 11-19 The 
basis for the strategy used here is suggested by the results of optimization work done by Miele and 
others. 20 The following is an interpretation of two important results from that work: 

1. Given atmospheric entry velocity and flight path angle, the optimal trajectory is a 2-arc 
trajectory; an entry arc flown at full positive lift, and an exit arc flown at full negative lift 

2. Given atmospheric entry velocity only, the optimal trajectory is a 1-arc trajectory flown at full 
negative lift 

Here, optimal is used in the sense of minimum peak heating rate, minimum structural loading, and 
minimum fuel usage for post-aeropass orbit insertion. Strictly speaking, the peak heating rate and 
structural loading are not mathematically minimized. However, optimization results have shown that the 
optimal trajectories discussed here are characterized by relatively small values of peak heating rate and 
structural loading. The optimal trajectories are constrained by equation (6), the apoapsis constraint, at 
aeropass exit Note that result 2 can be thought of as a special case of result 1 where the switch to Ml 
negative lift occurs at entry. 

Numerical generation of optimal trajectories in a nonreal-time environment is simple (once the 
form is known) in that in the case of result 1, all that is needed is to iterate on a switching time (i.e., the 
time to switch from full positive lift to full negative lift) such that the vehicle attains the target radius of 
apoapsis at aeropass exit In the case of result 2, all that is needed is to determine the shallowest flight 
path angle such that at full negative lift, the vehicle attains the target radius of apoapsis at aeropass exit. 

In a real-time guidance environment, result 2 does not have much practical use since the concern 
is with the trajectory from atmospheric entry to exit, and the entry flight path angle obviously cannot be 
influenced once at entry. Rather, result 2 is suggestive of a premission or preentry optimization (to 
determine entry flight path) that must be done in conjunction with a real-time guidance algorithm based 
on result 1. Thus, the guidance strategy discussed here is derived with result 1 as a starting point. 

Result 1 cannot be used directly because the magnitude of dispersions (atmospheric, commanded 
bank angle reversals, navigation errors, etc.) that would be encountered during the exit (negative lift) 
phase is unknown. If result 1 is used directly and, for example, the encountered atmospheric density is 
less than the atmospheric model predicted, there would not be enough lift capability to prevent a 
premature “skipout.” What is needed is some margin for uncertainties that the vehicle will encounter 
during the exit phase. This margin can be obtained by “tricking” the guidance into “thinking” the 
vehicle’s full negative-lift capability is less than what it really is. This is done by introducing a “bank 
margin,” a premission determined positive number, which is used during the entry phase modeling of the 
exit (i.e., negative lift) arc. This has the effect of causing the guidance logic to “think” it needs to start 
the exit phase earlier than it really needs to. Once the exit phase has been triggered, we simply solve for 
the bank angle required to hit the target apoapsis. 

It turns out that there is no need to rigorously “solve” for the switching time. All that is required 
is to generate one integrated trajectory per guidance cycle (during the entry phase) and monitor the 
results. The following algorithm is used as the entry phase guidance logic: 
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1. Given the current vehicle state and estimates for vehicle characteristics, numerically integrate 
the chosen equations of motion from to to cpt-Afjou using full lift up. Here, Afrou is the estimated time 
duration required for the vehicle to roll from full positive lift to full negative lift 

2. Given the simulated vehicle state from step 1, numerically integrate from fb+AJk>n» using a 
bank angle corresponding to full negative lift less the bank margin, until the simulated vehicle either 
becomes trapped in the atmosphere or exits the atmosphere. 

3. If the simulated vehicle becomes trapped in the atmosphere, then command the bank angle 
corresponding to full positive lift for the vehicle, i.e., when 


Sgn(C/J = 1 , <Pcmd = 0 , 


when 


sgn(Q)=-l , &md= 180. (14) 

4. If the simulated vehicle exits the atmosphere, then compute the radius of apoapsis. If the 
radius of apoapsis is smaller than the target radius of apoapsis, then command a full positive lift bank 
angle. If the radius of apoapsis is larger than the target, then the exit phase is initiated, and exit phase 
logic will be executed on subsequent guidance cycles. 

For the exit phase logic, simply use the basic apoapsis controller discussed above. 

A most pleasing feature to note here is that all of the gains and many of the mission-dependent 
parameters so common in other guidance algorithms have been replaced by a single mission-dependent 
parameter, i.e., the bank margin in step 2. Another elegant feature is that the time to transition from entry 
to exit phase is determined completely in real time and autonomously (no trigger or transition velocities 
are necessary). A computationally desirable feature is that in the entry phase, only one integrated trajec- 
tory need be generated and, furthermore, early in the aeropass these are short-lived trajectories as simu- 
lated vehicles become trapped in the atmosphere. The exit phase, which is inherently more computa- 
tionally intensive, does not start until the time-to-exit has been reduced somewhat in the entry phase. 


IV. REAL-TIME ORBITAL PLANE CONTROLLER 


The form of the plane controller is based on a first-order system where the state variable being 
controlled is the out-of-plane velocity component, that is, the vehicle’s velocity component along a unit 
vector perpendicular to the desired orbit plane. The only time this quantity is zero (during any finite time 
period) is when the vehicle’s orbit plane is identically the desired orbit plane. Thus, if the magnitude of 
out-of-plane velocity is driven to a small value, the same should automatically be accomplished with the 
plane error. Starting with the homogeneous equation: 

dy/dt+y/t=Q , (15) 

where y is the current out-of-plane velocity (measured positive along the unit vector antiparallel to the 
desired unit angular momentum vector), one assumes the out-of-plane lift component is the only signifi- 
cant out-of-plane force to obtain 
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La sin (fn-yA = 0 . 


(16) 


Solving for sin <f> and expanding L a one obtains 

sin l = (y) / rpcst C*£/(2 m) , (17) 

where fe *t is the estimated density at the current altitude. The sign of sin 0 gives the appropriate sign of 
the bank command, and the magnitude of sin 0 is a measure of the current plane error. The following 
defines the plane controller logic: 

When the magnitude of sin $ is larger than some positive number 5 max , set the sign of the 
commanded bank angle to that of sin <f> 


sgn(kmd) = sgn(sin 0) . (18) 

This commanded bank sign will be used in subsequent guidance cycles until sin ^ changes sign again 
and Isin <f > I becomes larger than S max . Inspection of equation (17) reveals that the magnitude of sin <f> 
becomes larger as dynamic pressure decreases. This results in increased sensitivity to plane error as 
atmospheric exit is approached and less sensitivity to plane error when dynamic pressure is large. This is 
precisely what is desired because large plane errors deep in the atmosphere are acceptable, i.e, there is 
plenty of controllability available during the ascent to atmospheric exit to correct for large plane errors. 
Note that the time constant, x, and S max are not independent parameters, but rather a doubling in x, for 
example, results in the same performance as does a doubling in S maT . Therefore, S max is arbitrarily set to 
one, and satisfactory guidance performance is achieved by adjusting x. Thus, as in the apoapsis con- 
troller case, the fortuitous situation occurs wherein only one critical parameter value needs to be 
(semi-)rigorously determined prior to a mission. An important point is that large values of x result in 
fewer roll reversals and large plane error at exit, while conversely, small values of x result in more roll 
reversals and small plane error at exit. Thus, it is a simple matter to tradeoff roll reversals for plane error 
to obtain satisfactory performance. 


V. NUMERICAL TESTING 


A version of the guidance logic discussed above has been coded into a FORTRAN subroutine, 
GREGOAER, in as general a manner as possible to enable use with different planets, vehicles, and entry 
conditions. A three degree-of-freedom (DOF) aerobraking guidance algorithm test-bed program was 
used to test both GREGOAER and a version of what was to be flown as AFE guidance. The test mission 
is essentially the AFE mission (geosynchronous Earth orbit (GEO) to low-Earth orbit (LEO) transfer 
with L/D ~ 0.3) with a nominal vehicle mass of 1,860 kg, nominal entry speed of 10,308 m/s, and nomi- 
nal entry flight path angle of -4.45°. No additional effort was made to optimize the AFE guidance per- 
formance over and above what was done during its extensive development Likewise, no effort was 
made to optimize the guidance performance of GREGOAER, except for the newly developed lateral 
guidance logic. A 45° bank margin was chosen for exit arc modeling in the GREGOAER entry phase. 
This choice was completely arbitrary. Vehicle mass dispersions (±13.6 kg, 1 o), entry flight path disper- 
sions (±0.01°, la) and atmospheric dispersions (generated using the GRAM atmosphere 21 ) were 
modeled for 100 test trajectories. 
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VI. NUMERICAL RESULTS 


Time histories of several important variables are shown in figures 1 through 7 for the nominal 
aeropass trajectory guided by GREGOAER. The high-frequency components seen in the relative velo- 
city, heat rate, and dynamic pressure plots are due to the perturbing winds generated by GRAM. The 
test-bed program provides a realistic modeling of the finite response time required for the actual bank 
angle to track the commanded bank angle. This is seen in the bank angle plot, figure 6. To avoid losing 




0 50 100 150 200 250 300 350 400 450 500 

Time From Reentry (s) 


Figure 2. Nominal trajectory: inertial velocity magnitude versus time. 
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Time From Reentry (s) 

Figure 3. Nominal trajectory: inertial flight path angle versus time. 
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Figure 4. Nominal trajectory: total heat rate versus time. 
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Time From Reentry (s) 

Figure 7. Nominal trajectory: bank angle versus time. 

plane controllability, the commanded bank angle magnitude is limited between 15° and 165°. This limit- 
ing typically occurs in the latter part of the exit phase and is seen in figure 6 after about 325 s into the 
aeropass. By limiting the bank angle, only about 3 percent of inplane lift capability is sacrificed, while a 
relatively large out-of-plane lift component (0.26 La) is provided. Roll reversals at about 145, 210, and 
310 s are seen in the bank angle plot The effect of these is seen in figure 7. As demonstrated in the out- 
of-plane velocity plot, the plane controller becomes more sensitive (more likely to command a roll 
reversal) to plane error as dynamic pressure, hence controllability, decreases. 

Statistics of important trajectory characteristics for the 100 trajectories using the AFE guidance 
are given in table 1 , while the results using GREGOAER are given in table 2. 


Table 1. Results using AFE guidance algorithm. 



Mean 

Min 

Max 

Minimum altitude, m 

75,138 

74,368 

75,808 

Peak G-load, g 

2.04 

1.85 

2.25 

Peak dyn pr, N/m 2 

1,800 

1,652 

2,025 

Peak ht rt, W/cm 2 

51.0 

47.6 

54.1 

Heat load, J/cm 2 

8,152 

7,768 

8,472 

Exit apogee, km 

328.0 

289.3 

344.1 

Exit perigee, km 

12.8 

-11.7 

42.0 

Exit wedge, degree (°) 

9.76E-2 

7.23D-2 

0.14 

Delta V, m/s 

106.8 

97.5 

118.3 

Roll Reversals 

2.8 

2 

4 
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Table 2. Results using GREGOAER guidance algorithm. 



Mean 

Min 

Max 

Minimum altitude, m 

76,133 

75,066 

77,178 

Peak G-load, g 

1.80 

1.57 

2.09 

Peak dyn pr, N/m 2 

1,594 

1,379 

1,853 

Peak ht rt, W/cm 2 

49.6 

46.0 

53.1 

Heat load, J/cm 2 

8,816 

8,382 

9,380 

Exit apogee, km 

342.2 

319.0 

359.4 

Exit perigee, km 

52.9 

25.8 

77.4 

Exit wedge, degree (°) 

3.7E-2 

1.1E-2 

6.4E-2 

Delta V, m/s 

90.1 

80.7 

98.8 

Roll Reversals 

3.9 

3 

5 


Note that the minimum altitudes experienced using GREGOAER are higher than those of AFE 
guidance resulting in lower loads and lower peak heat rates. GREGOAER total heat loads are higher 
because the times spent in the atmosphere are typically higher (by about 100 s). Neither AFE guidance 
nor GREGOAER had any problem hitting the target apogee of 340 km, although GREGOAER had a 
smaller spread between the maximum and minimum values which is a desirable feature. GREGOAER 
had higher exit perigees resulting in smaller post-aeropass delta V. GREGOAER did a better job mini- 
mizing the exit wedge angle (angle between the desired and actual planes), at a cost of one more roll 
reversal on average. A plane controller value of r= 45 s was used here. The average and worst-case cir- 
cularization (at 296 km) delta-V requirements for GREGOAER are significantly smaller due mainly to 
the higher exit perigees. 


VII. CONCLUSIONS 


With a bare minimum of development effort, the guidance algorithm discussed here is producing 
desirable results. Because it was based on sound physical principles and on optimality results, there is 
reason to believe it can provide superior and robust guidance performance for a wide variety of aero- 
braking missions with a minimum of laborious and costly premission guidance performance analyses. 
Current and future work in the area of onboard density profile estimators would also significantly 
improve guidance performance by reducing the bank margin required in the entry phase logic. 
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